home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / cmln1086.arc / HOLOGRAM.PAS < prev   
Pascal/Delphi Source File  |  1986-07-29  |  16KB  |  457 lines
  1.  
  2.  
  3.  
  4. {          PROGRAM:              ALIAS.PAS
  5.  
  6.            AUTHOR:               Brian Hayes
  7.  
  8.                                  Copyright 1986 Brian Hayes
  9.                                  May be freely copied for noncommercial use.
  10.  
  11.            VERSION:              1.0
  12.  
  13.            DATE BEGUN:           June 28, 1986
  14.  
  15.            LAST REVISION:        June 30, 1986
  16.  
  17.  
  18.            FOR COMPILATION BY:   Turbo Pascal v. 3.0 (8087 recommended)
  19.  
  20.            "INCLUDE" FILES:      None
  21. }
  22.  
  23.  
  24.  
  25. {          DESCRIPTION:
  26.  
  27.             Accepts a polynomial in two unknowns (X and Y) and plots
  28.             a contour map of Z as a function of X and Y over a specified
  29.             range of values. The contour map is printed on a dot-matrix
  30.             printer. The function is "rasterized" by evaluating it at
  31.             every point in the map and turning the corresponding pixel
  32.             on only if the resulting value, modulo the contour interval,
  33.             falls in a specified range. If the surface is steep enough
  34.             that there are fewer than two dots per contour interval,
  35.             aliasing results.
  36. }
  37.  
  38. {          COMPILER DIRECTIVES:     }
  39.  
  40. {$B+}   {B+ assigns StdIn/StdOut to CON, B- to TRM; default +}
  41. {$C-}   {C+ allows ^C and ^S during Read/ReadLn; default +}
  42. {$I+}   {I+ enables automatic I/O error checking; default +}
  43. {$R-}   {R+ enables run-time checking of index bounds; default -}
  44. {$V-}   {V+ requires string parameters to match declared length; default +}
  45. {$U-}   {U+ allows ^C interrupt at any time; default -}
  46. { $G0}  {G>0 allocates input buffer for I/O redirection; default 0}
  47. { $P0}  {P>0 allocates output buffer for I/O redirection; default 0}
  48. {$D+}   {D+ unbuffers I/O for devices; default +}
  49. {$F16}  {Fn sets maximum number of files open simultaneously; default 16}
  50. {$K-}   {K+ enables checking for stack-heap collision; default +}
  51.  
  52.  
  53.  
  54. program HOLOGRAM;
  55.  
  56. type
  57.   TermPtr = ^TermType;
  58.  
  59.   TermType = record           { Stores data on each term of the polynomial }
  60.     Coef     : real;
  61.     Xpower   : integer;
  62.     Ypower   : integer;
  63.     NextTerm : TermPtr;
  64.   end;
  65.  
  66.   RegPack = record
  67.     AL, AH, BL, BH, CL, CH, DL, DH : byte
  68.   end;
  69.  
  70.   InputStrType = string[80];
  71.  
  72. var
  73.   FirstTerm : TermPtr;          { Start of chain of terms  }
  74.   Regs      : RegPack;          { For MsDos call           }
  75.   MinXY     : real;             { Upper left corner        }
  76.   MaxXY     : real;             { Lower right corner       }
  77.   MapSize   : real;             { In inches; always square }
  78.   ContIntv  : real;             { Contour interval         }
  79.   ContWidth : real;             { Fraction of black/white  }
  80.   Istr      : InputStrType;     { Holds input equation     }
  81.   Answer    : char;             { For main program query   }
  82.  
  83.  
  84.  
  85.     { For a description of the following integer-power algorithm and
  86.       a historical account of its development see Dennis E. Hamilton,
  87.       "Fast Integer Powers for Pascal," Dr. Dobbs Journal, February,
  88.       1986, page 36. }
  89.  
  90.  
  91. function Power(N : real; E : integer): real;
  92.   var
  93.     TmpN : real;
  94.     TmpE : integer;
  95.   begin
  96.     case E of
  97.       0 :  Power := 1.0;   { NOTE : 0^0 error ignored }
  98.       1 :  Power := N;
  99.       2 :  Power := Sqr(N);
  100.       else begin
  101.              TmpE := Abs(E);
  102.              while not Odd(TmpE) do
  103.                begin
  104.                  N := Sqr(N); TmpE := TmpE shr 1;
  105.                end;
  106.              TmpN := N;
  107.              while TmpE > 1 do
  108.                begin
  109.                  repeat
  110.                    N := Sqr(N); TmpE := TmpE shr 1;
  111.                  until Odd(TmpE);
  112.                  TmpN := TmpN * N;
  113.                end;
  114.              if E < 0 then Power := 1.0/TmpN else Power := TmpN;
  115.            end;
  116.     end;
  117.   end;
  118.  
  119. function EvalPoly(X, Y : real) : real;
  120.  
  121.   { Evaluates the equation accepted by GetPoly at a single point. The
  122.     terms of the equation are stored in a linked list, which is traced
  123.     by EvalPoly. For each term the coefficient and the exponents stored
  124.     in the list node are applied, and the result is added to the running
  125.     total of all terms evaluated so far. (Negative terms are handled by
  126.     negating the coefficient in GetPoly.) When there are no more terms,
  127.     the function returns the total value.                                }
  128.  
  129.   var
  130.     Next  : TermPtr;
  131.     Value : real;
  132.   begin
  133.     Next := FirstTerm;
  134.     Value := 0.0;
  135.     while Next <> nil do
  136.       begin
  137.         if Next^.Coef <> 0.0 then Value :=
  138.           Value + (Next^.Coef * Power(X,Next^.Xpower) * Power(Y,Next^.Ypower));
  139.         Next := Next^.NextTerm;
  140.       end;
  141.     EvalPoly := Value;
  142.   end;
  143.  
  144. procedure GetPoly;
  145.  
  146.   { GetPoly accepts an equation from the user and parses it into terms;
  147.     each term consists of a coefficient and two exponents (for the X and
  148.     Y factors). If X or Y is absent from the term, an exponent of zero is
  149.     stored. If both X and Y are missing (i.e., the term is a constant),
  150.     two zero exponents are stored. Note that any nonzero value raised to
  151.     the zero power is 1.                                                   }
  152.  
  153.   { The format for an equation can be defined by the following example:
  154.                     Z = 2X^3Y^4 - 5XY + 6Y^7 - 8
  155.              Each term can have one coefficient, which must be the
  156.              first item in the term and must be recognizable by Turbo
  157.              Pascal as a real number. If the coefficient is omitted,
  158.              it is assumed to be 1. An exponents is an integer that
  159.              follows an X or a Y and a caret; missing exponents are
  160.              set to zero. Terms are separated by plus or minus signs.    }
  161.  
  162.   { WARNING: This routine is not robust. It is a quick-and-dirty
  163.              equation parser, which handles only simple polynomials
  164.              in X and Y. Useful extensions would include the ability
  165.              to recognize parenthesized expressions, the division
  166.              operator, fractional exponents and trig functions. It
  167.              would also be convenient to enter equations in polar
  168.              coordinates. Finally, the error handling included here
  169.              is very crude.                                           }
  170.  
  171.   type
  172.     PolyStrType = string[82];
  173.     ErrorCode = (BadCoef, BadExp);
  174.   var
  175.     Pstr : PolyStrType;       { Holds equation string }
  176.     This : TermPtr;           { Current term pointer  }
  177.     Old  : TermPtr;           { Needed to forge links }
  178.     I    : 0..83;             { Index into Pstr    }
  179.     ParseMore : boolean;      { Flag to exit loop  }
  180.     ErrorFlag : boolean;      { Flag to mark error }
  181.  
  182.   procedure Error(Err: ErrorCode);
  183.     begin
  184.       ParseMore := False;                             { Exit GetPoly          }
  185.       if not ErrorFlag then                           { Report 1st error only }
  186.         begin
  187.           ErrorFlag := True;
  188.           case Err of
  189.             BadCoef : WriteLn('Bad coefficient.');
  190.             BadExp  : WriteLn('Bad exponent.');
  191.           end;
  192.         end;
  193.     end;
  194.  
  195.   function GetChar : char;            { Get and remove next char from string }
  196.     begin
  197.       I := Succ(I);                   { Must not be called at end of string }
  198.       GetChar := Pstr[I];
  199.     end;
  200.  
  201.   function NextChar : char;           { Peek at next but don't remove }
  202.     begin
  203.       NextChar := Pstr[Succ(I)];
  204.     end;
  205.  
  206.   procedure GetCoef;
  207.     var
  208.       CoefSign : char;
  209.       CoefStr  : InputStrType;
  210.       CoefVal  : real;
  211.       ConvCode : integer;
  212.     begin
  213.       if not (NextChar in ['+','-']) then Error(BadCoef);    { Must have sign }
  214.       while (NextChar in ['+','-']) do CoefSign := GetChar;  { Use last sign  }
  215.       CoefStr := '';
  216.       while (NextChar in ['0'..'9','.']) do CoefStr := CoefStr + GetChar;
  217.       if (CoefStr = '') then This^.Coef := 1.0
  218.         else
  219.           begin
  220.             Val(CoefStr, CoefVal, ConvCode);
  221.             if ConvCode = 0 then This^.Coef := CoefVal else Error(BadCoef);
  222.           end;
  223.       if CoefSign = '-' then This^.Coef := -This^.Coef;
  224.       if (NextChar = #00) then ParseMore := False;        { End of string }
  225.     end;
  226.  
  227.   procedure GetX;
  228.     var
  229.       Dump     : char;
  230.       ExpSign  : char;
  231.       ExpStr   : InputStrType;
  232.       ExpVal   : integer;
  233.       ConvCode : integer;
  234.     begin
  235.       if (NextChar in ['X','x']) then
  236.         begin
  237.           Dump := GetChar;
  238.           if (NextChar = '^') then
  239.             begin
  240.               Dump := GetChar;
  241.               while (NextChar in ['+','-']) do ExpSign := GetChar;
  242.               ExpStr := '';
  243.               while (NextChar in ['0'..'9']) do ExpStr := ExpStr + GetChar;
  244.               if ExpStr = '' then Error(BadExp);
  245.               Val(ExpStr, ExpVal, ConvCode);
  246.               if ConvCode = 0 then This^.Xpower := ExpVal else Error(BadExp);
  247.               if ExpSign = '-' then This^.Xpower := -This^.Xpower;
  248.             end
  249.           else This^.Xpower := 1
  250.         end;
  251.       if (NextChar = #00) then ParseMore := False;
  252.     end;
  253.  
  254.   procedure GetY;
  255.     var
  256.       Dump     : char;
  257.       ExpSign  : char;
  258.       ExpStr   : InputStrType;
  259.       ExpVal   : integer;
  260.       ConvCode : integer;
  261.     begin
  262.       if (NextChar in ['Y','y']) then
  263.         begin
  264.           Dump := GetChar;
  265.           if (NextChar = '^') then
  266.             begin
  267.               Dump := GetChar;
  268.               while (NextChar in ['+','-']) do ExpSign := GetChar;
  269.               ExpStr := '';
  270.               while (NextChar in ['0'..'9']) do ExpStr := ExpStr + GetChar;
  271.               if ExpStr = '' then Error(BadExp);
  272.               Val(ExpStr, ExpVal, ConvCode);
  273.               if ConvCode = 0 then This^.Ypower := ExpVal else Error(BadExp);
  274.               if ExpSign = '-' then This^.Ypower := -This^.Ypower;
  275.             end
  276.           else This^.Ypower := 1
  277.         end;
  278.       if (NextChar = #00) then ParseMore := False;
  279.     end;
  280.  
  281.   procedure NewTerm;
  282.     begin
  283.       Old := This;               { Save pointer so we can add link }
  284.       New(This);
  285.       Old^.NextTerm := This;     { Link to the new record          }
  286.       This^.Coef   := 0.0;
  287.       This^.Xpower := 0;
  288.       This^.Ypower := 0;
  289.       This^.NextTerm := nil;     { Mark this as last term in chain }
  290.     end;
  291.  
  292.   begin   { GetPoly }
  293.     repeat
  294.       ErrorFlag := False;
  295.       WriteLn('Enter a polynomial in X and Y (e.g., Z = 2X^2 - Y + 5):');
  296.       WriteLn;
  297.       Write('Z = ');
  298.       ReadLn(Istr);               { Get equation from terminal }
  299.       WriteLn; WriteLn;
  300.       Pstr := '+' + Istr + #00;   { Precede first term with plus sign }
  301.       for I := 1 to Length(Pstr) do                   { Remove spaces }
  302.        if (Pstr[I] = #32) then Delete(Pstr, I, 1);
  303.       Mark(This);
  304.       NewTerm;
  305.       FirstTerm := This;          { Set pointer to start of chain }
  306.       I := 0;
  307.       ParseMore := True;
  308.       repeat
  309.         if ParseMore then GetCoef;
  310.         if ParseMore then GetX;
  311.         if ParseMore then GetY;
  312.         if ParseMore then NewTerm;
  313.       until not ParseMore;
  314.     until ErrorFlag = False;
  315.   end;
  316.  
  317. procedure GetParams;
  318.   begin
  319.     repeat
  320.       Write('Range of X and Y to be plotted: From ');
  321.       Read(MinXY);
  322.       Write(' to ');
  323.       ReadLn(MaxXY);
  324.       WriteLn;
  325.     until (MaxXY <> MinXY);
  326.     repeat
  327.       Write('Map size (inches): ');
  328.       ReadLn(MapSize);
  329.       WriteLn;
  330.     until (MapSize > 0.25);
  331.     repeat
  332.       Write('Contour interval: ');
  333.       ReadLn(ContIntv);
  334.       WriteLn;
  335.     until (ContIntv > 0);
  336.     repeat
  337.       Write('Contour width (between 0.0 and 1.0): ');
  338.       ReadLn(ContWidth);
  339.       WriteLn;
  340.     until ((ContWidth > 0.0) and (ContWidth < 1.0));
  341.   end;
  342.  
  343. function GetInput: boolean;
  344.   var
  345.     Yes : char;
  346.   begin
  347.     GetInput := False;
  348.     ClrScr;
  349.     WriteLn('CONTOUR MAPPING');
  350.     WriteLn; WriteLn;
  351.     GetPoly;
  352.     GetParams;
  353.     Write('Proceed to plot map? (Y/N) ');
  354.     repeat
  355.       Read(Kbd, Yes);
  356.       if UpCase(Yes) = 'Y' then GetInput := True;
  357.     until UpCase(Yes) in ['Y','N'];
  358.   end;
  359.  
  360. procedure PrintByte(Pattern : byte);
  361.  
  362.   { The Pascal WRITE procedure does not accept a parameter of type byte.
  363.     This call on MsDos outputs a byte directly to the printer. Changes
  364.     will be needed for other operating systems.                           }
  365.  
  366.   begin
  367.     Regs.AH := 5;             { Printer output service }
  368.     Regs.DL := Pattern;       { Byte to be printed     }
  369.     MsDos(Regs);
  370.   end;
  371.  
  372.  
  373.  
  374.   { NOTE: The following procedures include escape codes for the
  375.           graphics modes of the Epson LQ-1500 printer. Other
  376.           Epson printers use similar (if not identical codes).     }
  377.  
  378.  
  379. procedure PrintBorder(N1, N2, Pattern : byte);
  380.  
  381.   { Prints a full line of a single graphics pattern. Used to produce
  382.     a single row of dots top and bottom as a border for the map.      }
  383.  
  384.   var
  385.     Row : integer;
  386.   begin
  387.     Write(Lst,#27'*'#0);                { Graphics mode zero       }
  388.     PrintByte(N1); PrintByte(N2);       { Number of bytes to print }
  389.     for Row := 0 to (Integer(N2*256+N1)) do PrintByte(Pattern);
  390.     WriteLn(Lst);
  391.   end;
  392.  
  393. procedure PrintInfo;
  394.   begin
  395.     Write(Lst, 'Z = ');
  396.     WriteLn(Lst, Istr);
  397.     WriteLn(Lst);
  398.     WriteLn(Lst, 'From:     ', MinXY:10:2);
  399.     WriteLn(Lst, 'To:       ', MaxXY:10:2);
  400.     WriteLn(Lst, 'Interval: ', ContIntv:10:2);
  401.     WriteLn(Lst, 'Width:    ', ContWidth:10:2);
  402.     WriteLn(Lst);
  403.     WriteLn(Lst);
  404.   end;
  405.  
  406. procedure PlotMap;
  407.   var
  408.     I, J, K  : integer;    { Loop indices                       }
  409.     X, Y     : real;       { Map coordinates                    }
  410.     Delta    : real;       { Change in X or Y per dot interval  }
  411.     Dots     : integer;    { Number of dots per line            }
  412.     N1, N2   : byte;       { Number of dots in Epson format     }
  413.     ThisByte : byte;       { Holds pattern as points calculated }
  414.   begin
  415.     Dots := (Trunc(MapSize * 60) and $FFF8);           { (div 8) cols & rows  }
  416.     N1 := Byte(Lo(Dots + 2));                          { Cols per row, in...  }
  417.     N2 := Byte(Hi(Dots + 2));                          { ... Epson format     }
  418.     Delta := (MaxXY - MinXY) / Dots;                   { Z(X,Y) sample points }
  419.     PrintInfo;                                         { Identify the plot    }
  420.     WriteLn(Lst,#27'3'#24);                            { Set 24/180 linespace }
  421.     PrintBorder(N1, N2, 1);                            { Top border           }
  422.     for K := 0 to ((Dots div 8) - 1) do                { Each K one line      }
  423.       begin
  424.         Write(Lst,#27'*'#0);                           { 60 dpi graphics      }
  425.         PrintByte(N1); PrintByte(N2);                  { Number of bytes      }
  426.         PrintByte(255);                                { Left border          }
  427.         for J := 0 to (Dots - 1) do                    { Each J one column    }
  428.           begin
  429.             ThisByte := 0;                             { Initially all blank  }
  430.             for I := 0 to 7 do                         { Each I one dot       }
  431.               begin
  432.                 X := MinXY + (J * Delta);                              {      }
  433.                 Y := MinXY + ((K * 8 + I) * Delta);                    { SEE  }
  434.                 if (Frac(Abs((EvalPoly(X,Y))/ContIntv)) < ContWidth)   { NOTE }
  435.                   then ThisByte := ThisByte or 128 shr I;              {      }
  436.               end;
  437.             PrintByte(ThisByte);                       { Print one column     }
  438.           end;
  439.         PrintByte(255);                                { Right border         }
  440.         WriteLn(Lst);                                  { CR/LF to printer     }
  441.         if Keypressed then begin Read(Kbd); exit end;  { Allow user break     }
  442.       end;
  443.     PrintBorder(N1, N2, 128);                          { Bottom border        }
  444.     WriteLn(Lst,#27'2');                               { Reset line spacing   }
  445.   end;
  446.  
  447.  
  448.  
  449. begin                { Main Program }
  450.   repeat
  451.     if GetInput then PlotMap;
  452.     WriteLn; WriteLn; WriteLn;
  453.     Write('<P>lot another.    <Q>uit.');
  454.     repeat Read(Kbd, Answer) until UpCase(Answer) in ['P','Q'];
  455.   until UpCase(Answer) = 'Q';
  456. end.
  457.